########################################################################
# PLS 206 — Week 3 live coding, part 2: generalized linear models
# Grey Monroe, Fall 2026
########################################################################

# install.packages("ggplot2")
library(ggplot2)

trial <- read.csv("https://greymonroe.github.io/PLS_206/data/wheat_trial.csv")
head(trial)

########################################################################
# 1. A straight line fails for yes/no
########################################################################

# rust: 1 = stripe rust on the plot, 0 = none

m_lin <- lm(rust ~ nitrogen + variety + water, data = trial)
range(fitted(m_lin))                          # fitted "probabilities" for the real plots
sum(fitted(m_lin) < 0 | fitted(m_lin) > 1)   # how many are impossible

########################################################################
# 2. Probability, odds, log-odds
########################################################################

p <- c(0.2, 0.5, 0.8)
odds <- p / (1 - p)
cbind(p, odds, log_odds = log(odds))
qlogis(0.8)    # probability -> log-odds
plogis(1.39)   # log-odds -> probability

########################################################################
# 3. Logistic regression with a group
########################################################################

table(trial$rust, trial$water)
m_rw <- glm(rust ~ water, family = binomial, data = trial)
summary(m_rw)                 # log-odds scale

plogis(coef(m_rw)[1])         # P(rust), control
plogis(sum(coef(m_rw)))       # P(rust), drought

exp(coef(m_rw))               # odds ratio
exp(confint(m_rw))

########################################################################
# 4. Logistic regression with a slope
########################################################################

m_rn <- glm(rust ~ nitrogen, family = binomial, data = trial)
summary(m_rn)
exp(10 * coef(m_rn)["nitrogen"])   # odds ratio per 10 kg N

# predict(): log-odds unless you ask for the response
new <- data.frame(nitrogen = c(50, 100, 150))
predict(m_rn, newdata = new)
predict(m_rn, newdata = new, type = "response")

# the S-curve
grid <- data.frame(nitrogen = 0:200)
grid$p <- predict(m_rn, newdata = grid, type = "response")
ggplot(trial, aes(nitrogen, rust)) +
  geom_jitter(width = 0, height = 0.03) +
  geom_line(data = grid, aes(nitrogen, p), colour = "red") +
  theme_classic()

########################################################################
# 5. How often is it right?
########################################################################

pred <- as.numeric(fitted(m_rn) > 0.5)
table(predicted = pred, observed = trial$rust)
mean(pred == trial$rust)      # accuracy
mean(trial$rust == 0)         # always guess "no rust"

########################################################################
# 6. Several predictors and deviance
########################################################################

m_rall <- glm(rust ~ nitrogen + variety + water, family = binomial, data = trial)
summary(m_rall)
round(exp(coef(m_rall)), 2)   # odds ratios

anova(m_rn, m_rall, test = "Chisq")   # do variety and water add anything?

########################################################################
# 7. Counts: Poisson regression
########################################################################

hist(trial$aphids)
m_aph <- glm(aphids ~ nitrogen, family = poisson, data = trial)
summary(m_aph)
exp(10 * coef(m_aph)["nitrogen"])     # rate ratio per 10 kg N

# overdispersion: about 1 is fine, well above 1.5 is a problem
deviance(m_aph) / df.residual(m_aph)

########################################################################
# 8. Fix: negative binomial
########################################################################

library(MASS)
m_nb <- glm.nb(aphids ~ nitrogen, data = trial)
summary(m_nb)                         # same slope, honest SE

exp(10 * coef(m_nb)["nitrogen"])
exp(10 * confint(m_nb)["nitrogen", ])

########################################################################
# 9. Bonus: mixed effect models
########################################################################

# 112 plots on 12 farms (simulated). Cover crops were chosen per farm
# and in truth do nothing.

set.seed(181)
nj <- c(2, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 25)   # plots per farm
ids <- sprintf("farm%02d", 1:12)
cover_j <- setNames(rep(c("no", "yes"), 6)[sample(12)], ids)
u <- setNames(rnorm(12, 0, 0.8), ids)              # each farm: its own baseline
farm <- rep(ids, nj)
nitrogen <- round(runif(112, 0, 200))
yield <- round(3.5 + 0.012 * nitrogen + u[farm] + rnorm(112, 0, 1), 2)
farms <- data.frame(farm, cover = unname(cover_j[farm]), nitrogen, yield)

# install.packages("lme4")
library(lme4)
m_mixed <- lmer(yield ~ nitrogen + (1 | farm), data = farms)
VarCorr(m_mixed)        # how much farms differ vs plots within a farm
summary(m_mixed)

# a treatment chosen per farm: plots are not independent evidence
summary(lm(yield ~ nitrogen + cover, data = farms))
summary(lmer(yield ~ nitrogen + cover + (1 | farm), data = farms))